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A  non-equilibrium  molecular  dynamics  method  is  employed  to  study  the  mechanical  response  of  soda-lime 
glass  (a  material  commonly  used  in  transparent  armor  applications)  when  subjected  to  the  loading  con¬ 
ditions  associated  with  the  generation  and  propagation  of  planar  shock  waves.  Specific  attention  is  given  to 
the  identification  and  characterization  of  various  (inelastic-deformation  and  energy-dissipation)  molecular- 
level  phenomena  and  processes  taking  place  at,  or  in  the  vicinity  of,  the  shock  front.  The  results  obtained 
revealed  that  the  shock  loading  causes  a  2-4%  (shock  strength-dependent)  density  increase.  In  addition,  an 
increase  in  the  average  coordination  number  of  the  silicon  atoms  is  observed  along  with  the  creation  of 
smaller  Si-O  rings.  These  processes  are  associated  with  substantial  energy  absorption  and  dissipation  and 
are  believed  to  greatly  influence  the  blast/ballistic  impact  mitigation  potential  of  soda-lime  glass.  The 
present  work  was  also  aimed  at  the  determination  of  the  shock  Hugoniot  (i.e.,  a  set  of  axial  stress  vs.  density/ 
specific-volume  vs.  internal  energy  vs.  particle  velocity  vs.  temperature)  material  states  obtained  in  soda- 
lime  glass  after  the  passage  of  a  shock  wave  of  a  given  strength  (as  quantified  by  the  shock  speed).  The 
availability  of  a  shock  Hugoniot  is  critical  for  construction  of  a  high  deformation-rate,  large-strain,  high 
pressure  material  model  which  can  be  used  within  a  continuum-level  computational  analysis  to  capture  the 
response  of  a  soda-lime  glass  based  laminated  transparent  armor  structure  (e.g.,  a  military  vehicle  wind¬ 
shield,  door  window,  etc.)  to  blast/ballistic  impact  loading. 
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1.  Introduction 

Transparent  blast/ballistic-impact  resistant  vehicle  structures 
(e.g.,  windshields,  door  windows,  viewports,  etc.),  employed 
today,  utilize  a  variety  of  non-glass  transparent  materials  such 
as  transparent  crystalline  ceramics  (e.g.,  aluminum-oxinitride 
spinel,  AlON,  sapphire  (Ref  1)  and  new  transparent  polymer 
materials  (e.g.,  transparent  nylon  (Ref  2)).  Despite  the  well- 
established  benefits  offered  by  these  materials,  ballistic  glass 
remains  an  important  constituent  material  in  a  majority  of 
transparent  impact  resistant  structures  used  today.  The  main 
reason  for  this  are  as  follows:  (a)  glass-structure  fabrication 
technologies  enable  the  production  of  curved,  large  surface- 
area,  transparent  structures  with  thickness  approaching  several 
inches;  (b)  relatively  low  material  and  manufacturing  costs;  and 
(c)  compositional  modifications,  chemical  strengthening,  and 
controlled  crystallization  have  demonstrated  the  capability  to 
significantly  improve  the  shock/ballistic  impact  survivability  of 
glass  (Ref  2). 
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Typically,  the  development  of  new  glass-based  transparent 
impact  resistant  structures,  aimed  at  reducing  the  vulnerability 
of  protected  vehicle  occupants  and  on-board  instrumentation  to 
various  blast/ballistic  threats,  requires  extensive  prototyping 
and  laboratory/field  experimental  testing.  These  experimental 
efforts  are  critical  for  ensuring  the  effectiveness  and  reliability 
of  the  transparent  impact  resistant  structures. 

However,  these  efforts  are  generally  expensive,  time-con¬ 
suming,  and  involve  destructive  test  procedures.  While  the  role 
of  prototyping/testing  programs  remains  critical,  they  are 
increasingly  being  complemented  by  the  corresponding  com¬ 
putation-based  modeling  and  simulation  efforts.  It  is  well- 
established  (Ref  3-5),  that  the  effectiveness  and  reliability  of 
the  computation-based  modeling  and  simulation  approaches  is 
greatly  affected  by  the  fidelity  of  the  associated  material 
models.  It  is  critical  that  these  models  realistically  describe 
deformation/fracture  response  of  the  subject  material  (ballistic 
glass,  in  the  present  case)  under  high-rate,  large-strain,  high- 
pressure  loading  conditions  encountered  during  blast/ballistic 
impact.  Therefore,  one  of  the  main  objectives  of  the  present 
work  is  to  further  advance  the  application  of  computational 
modeling/simulation-based  engineering  approaches  of  trans¬ 
parent  impact-resistant  structures  via  the  identification  and 
quantification  of  processes  and  phenomena  occurring  in  glass 
under  such  conditions.  Specifically,  phenomena  and  processes 
associated  with  the  formation  and  propagation  of  shock  waves 
within  soda-lime  glass  will  be  investigated. 

A  comprehensive  literature  review  carried  out  as  part  of  our 
prior  work  (Ref  6)  revealed  that  the  mechanical  behavior  of 
glass  is  modeled  predominantly  using  three  distinct  approaches: 
(a)  molecular-modeling  methods  (Ref  7-10);  (b)  continuum- 
material  approximations  (Ref  3-5,  11-14);  and  (c)  models 
based  on  explicit  crack  representation  (Ref  15,  16).  Since  a 
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detailed  overview  of  these  models  can  be  found  in  Ref  6,  they 
will  not  be  discussed  any  further  in  the  present  manuscript. 
However,  what  is  worth  mentioning  is  that  the  overview 
presented  in  Ref  6  clearly  established  that:  (a)  molecular-level 
models  are  critical  for  identifying  nanometer  length  scale 
phenomena  and  the  associated  micro  structure  evolution  pro¬ 
cesses;  (b)  continuum-level  models  are  the  only  ones  which 
enable  large-scale  computational  analysis  of  the  mechanical 
response  of  real  transparent-armor  protective  structures  to 
various  blast  and  ballistic  threats;  and  (c)  the  models  based  on 
explicit  crack  propagation  are  mainly  suitable  for  modeling 
local  deformation  and  fracture  response  of  simple-geometry  test 
coupons.  Since  the  main  objective  of  the  present  work  is  to  help 
identify  and  characterize  various  phenomena  and  processes 
accompanying  shock  generation  and  propagation  within  soda- 
lime  glass,  only  the  molecular-level  computational  models  and 
methods  will  be  employed. 

As  mentioned  above,  the  main  objective  of  the  present  work 
is  to  investigate  various  shock  wave-related  phenomena  using 
molecular-level  computational  methods. 

A  shock  wave  (or  simply  a  shock)  is  a  wave  which 
propagates  through  a  medium  at  a  speed  higher  than  the  sound 
speed  and  its  passage  causes  an  abrupt  and  discontinuous 
change  in  the  material  state  variables  (e.g.,  pressure,  internal 
energy,  density,  temperature,  and  particle  velocity).  The 
magnitude  of  the  state-variable  changes  and  the  shock  speed 
increase  with  the  strength  of  the  shock.  While  acoustic  waves 
give  rise  to  isentropic  changes  in  the  material  state  variables, 
passage  of  a  shock  is  typically  associated  with  irreversible 
(entropy-increasing)  changes  in  the  same  variables.  The  reason 
behind  this  difference  is  that  shock  involves  very  high  strain 
rates  that  bring  energy-dissipative  viscous  effects  into  promi¬ 
nence. 

A  review  of  the  literature  carried  out  as  a  part  of  the  present 
work  revealed  that  the  molecular-level  computational  methods 
were  first  employed  to  study  shock  waves  more  than  30  years 
ago  (Ref  17-20).  While  the  initial  studies  focused  on  shock 
phenomena  in  dense  fluids,  subsequent  work  also  included 
single-crystalline  solids.  The  key  findings  related  to  the 
generation  and  propagation  of  shocks  in  these  solids  can  be 
summarized  as  follows:  (a)  these  phenomena  are  inherently 
more  complex  in  solids  than  in  fluids,  because  solids,  in 
addition  to  the  lattice  parameter,  introduce  a  new  length  scale 
(i.e.,  the  size/spacing  of  defects)  which  tends  to  control  the 
nature/extent  of  inelastic-deformation  and  dissipative  processes 
at  the  shock  front;  (b)  shock  propagation  typically  results  in  the 
formation  of  steady  (time-invariant)  waves,  and  this  process  is 
facilitated  by  the  transverse  displacement  of  atoms  which 
produce  inelastic  deformation.  This  deformation  involves 
concerted  slippage  of  atoms  over  each  other  (and  not  viscous 
flow  as  in  the  case  of  shocks  in  fluids);  (c)  in  order  to  eliminate 
free-surface  effects,  molecular-level  modeling  of  shock  is 
typically  carried  out  using  computational  systems  with  periodic 
boundary  conditions  (at  least  in  directions  transverse  to  the 
shock-wave  propagation  direction).  To  attain  steady  wave 
conditions  of  the  shock,  computational  domains  sufficiently 
long  in  the  direction  of  shock-wave  propagation  have  to  be 
employed.  Furthermore,  lateral  dimensions  of  the  computa¬ 
tional  domain  have  to  be  also  sufficiently  large  to  avoid 
spurious  effects  associated  with  the  use  of  the  periodic 
boundary  conditions.  Consequently,  computational  domains 
involving  several  tens  to  hundred  thousands  of  atoms  have  to 
be  employed.  The  corresponding  computational  times  (con¬ 


trolled  by  the  shock-wave  transversal  time)  are  then  on  the 
order  of  5-20  ps.  This  limits  the  shock  thicknesses  to  around  10 
nm  and  the  corresponding  rise  times  to  ca.  0.5-1  ps.  Thus,  weak 
shock  waves  with  thicknesses  of  hundreds  of  nanometers  could 
not  be  analyzed  using  molecular-level  methods  (at  least  in  their 
steady-wave  regime). 

While  recognizing  the  aforementioned  aspects  and  potential 
limitations  of  molecular-level  modeling  of  shock,  a  preliminary 
computational  study  of  shock  generation/propagation  in  soda- 
lime  glass  is  carried  out.  There  are  two  main  objectives  of  the 
present  work: 

(a)  To  determine  the  shock  Hugoniot  (centered  on  the  initial 
stress-free  quiescent  state)  of  soda-lime  glass.  A  Hugoni¬ 
ot  is  the  locus  of  axial  stress  vs.  specific  volume  vs.  en¬ 
ergy  density  (vs.  particle  velocity  vs.  shock  speed) 
shocked  (“upstream”)  material  states.  The  Hugoniot  is 
often  used  in  the  derivation  of  the  continuum-level  mate¬ 
rial  models  (particularly  in  the  derivation  of  the  equation 
of  state)  which  is  used  in  the  computational  investiga¬ 
tions  of  the  response  of  structures  to  shock  loading.  In 
situations  in  which  one  is  interested  only  in  the  problem 
of  planar  shock  propagation/interaction  (in  the  presence 
of  uni-axial  strain  deformation  states),  a  complete  defini¬ 
tion  of  the  continuum-level  material  model  is  not  re¬ 
quired.  Instead,  the  knowledge  of  the  corresponding 
Hugoniot  relations  is  sufficient;  and 

(b)  The  Hugoniot  relations  mentioned  in  (a)  provide  a  glo¬ 
bal  statement  of  mass,  momentum,  and  energy  conserva¬ 
tion  accompanying  shock-induced  material  transition 
from  a  given  initial  (“downstream”)  equilibrium  state  to 
all  possible  final  (“upstream”)  equilibrium  states  for 
steady  planar  shock  waves  (of  different  strengths).  How¬ 
ever,  these  relations  provide  no  information  about  the 
structure  of  the  shock  front  or  the  nature  of  the  dissipa¬ 
tive  structural  rearrangement  mechanisms  that  lead  to  a 
steady  shock  wave.  Hence,  the  second  objective  of  the 
present  work  is  to  carry  out  a  detailed  examination  of 
the  downstream,  shock-front,  and  upstream  material 
states  (as  represented  by  the  local  stresses,  strains,  densi¬ 
ties/specific  volumes,  temperatures,  etc.)  and  molecular- 
level  morphology  to  identify  and  characterize  these  pro¬ 
cesses. 

The  organization  of  the  paper  is  as  follows:  A  brief 
description  of  the  molecular-level  micro  structure  of  glass 
including  its  random-network  representation  is  presented  in 
Sect.  2.  A  brief  overview  of  the  molecular-level  computational 
procedure  including  the  computational  cell  construction,  force 
field  identification,  computational  method(s)  selection,  shock- 
wave  generation,  and  the  problem  definition  are,  respectively, 
presented  in  Sects.  3. 1-3.5.  The  key  results  obtained  in  the 
present  work  are  presented  and  discussed  in  Sect.  4.  A 
summary  of  the  work  carried  out  and  the  key  results/ 
conclusions  are  given  in  Sect.  5. 


2.  Molecular-Level  Microstructure  of  Glass 

One  of  the  basic  properties  of  all  types  of  glass  is  their  lack 
of  long-range  atomic  order  which  classifies  them  as  amorphous 
materials.  For  instance,  the  atomic  arrangement  in  pure  silicate 
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glass  (i.e.,  fused  silica)  is  highly  random  relative  to  the 
chemically  equivalent  (crystalline)  quartz.  To  describe  the 
random  atomic  arrangement  within  glass  the  so-called  “random 
network  model”  (Ref  21)  is  typically  employed.  Such  a  model 
represents  amorphous  materials  as  a  three-dimensional  linked 
network  of  polyhedra  with  central  cations  of  various  coordi¬ 
nation  depending  on  the  character  of  the  atomic  constituents.  In 
the  case  of  silicate-based  glasses  like  fused  silica  and  soda-lime 
glass,  the  dominant  cation  is  silicon.  In  this  case,  each  silicon  is 
surrounded  by  four  oxygen  anions  and  forms  a  Si044- 
tetrahedron,  whereby  each  oxygen  is  bonded  to  (or  bridges) 
two  silicon  atoms.  Other  polyhedra  may  exist  in  the  network 
depending  on  the  glass  additives.  Since  silicon  has  a  tendency 
to  form  a  continuous  network  with  (bridging)  oxygen  atoms, 
Si02  is  commonly  referred  to  as  a  “network  former”.  Other 
potential  network  formers  in  glass  are  boron  and  germanium 
oxides. 

Numerous  oxides  and  other  additives  are  used  to  modify  the 
basic  silica  tetrahedra  network  of  silicate-based  glasses  to  tailor 
their  properties  to  specific  applications.  When  alkali  (or  alkaline 
earth)  oxides  are  added  to  a  pure  silicate-based  glass,  in  order  to 
accommodate  the  excess  of  oxygen  anions  which  are  present 
due  to  the  oxide  dissociation,  the  continuity  of  the  silica 
tetrahedra  network  becomes  disrupted.  The  resulting  glass 
structure  contains  additional  non-bridging  (connected  to  only 
one  silicon  atom)  oxygen  atoms.  Charge  transfer  from  the  alkali 
earth  metal  atoms  converts  the  non-bridging  oxygen  atoms  into 
singly  charged  anions.  The  metallic  cations  formed  in  this 
process  tend  to  hover  around  the  non-bridging  oxygen  ions  for 
local  charge  neutrality.  Since  alkali  (or  alkaline  earth)-based 
oxides  cause  a  disruption  in  the  continuous  glass  network,  they 
are  typically  referred  to  as  “network  modifiers”.  In  the  specific 
case  of  soda-lime  glass,  the  14  wt.%Na20  and  9  wt.%CaO 
additions  to  the  silica  both  act  as  network  modifiers. 

In  contrast  to  network  modifying  oxides,  more  covalently 
bonded  oxides  like  B203  donate  metallic  cations  which  become 
incorporated  into  the  glass  network  without  significant  com¬ 
promise  in  the  network  connectivity.  In  other  words,  the 
addition  of  these  oxides  does  not  lead  to  the  formation  of  non¬ 
bridging  oxygen  atoms.  Due  to  the  aforementioned  behavior  of 
these  glass-oxides  additives  they  are  typically  referred  to  as 
network  formers. 

For  the  purpose  of  quantitative  description  of  the  micro¬ 
structure  within  the  random  network  model,  a  number  of 
network  parameters  have  been  introduced.  For  instance,  a  so- 
called  R  value  (defined  as  the  average  number  of  oxygen  ions 
per  network  forming  ion)  is  used  to  describe  the  overall 
connectivity  of  a  given  network.  In  the  case  of  fused  silica,  in 
which  there  are  two  (bridging)  oxygen  anions  for  every 
network  forming  silicon  cation,  the  R  value  is  2.0.  In  the  case  of 
soda-lime  glass,  the  introduction  of  non-bridging  oxygen  ions 
causes  the  R  value  to  increase  to  ca.  2.41.  This  finding  suggests 
that  a  higher  R  value  is  associated  with  a  more  open  (less 
connected)  weaker  amorphous  glass  network.  As  far  as  the 
network  formers  are  concerned  their  effects  on  the  glass 
micro  structure  and  the  R  value  is  greatly  affected  by  the  cation 
coordination  number  and  the  strength  of  its  bond  with  oxygen 
as  well  as  to  the  concentration  of  the  additive.  In  addition  to  the 
R  value,  the  so-called  “A”  and  “7”  parameters  are  often  used 
to  further  quantify  the  glass  network  structure.  These  two 
parameters  are,  respectively,  defined  as  the  average  number  of 
non-bridging  (connected  to  only  a  single  network  forming 
cation)  and  bridging  (connected  to  two  network  forming 


cations)  oxygen  atoms  per  network  polyhedron.  Since  pure 
silicate-based  glasses  contain  only  bridging  oxygen  atoms  (and 
are  based  on  silica  tetrahedra),  the  X  and  Y  parameters  take  on 
values  of  0.0  and  4.0,  respectively.  On  the  other  hand,  since 
soda-lime  glass  contains  additional  non-bridging  oxygen  atoms 
(while  its  network  is  still  based  on  silica  tetrahedral),  X  takes  on 
a  non-zero  value  (ca.  0.81)  while  Y  drops  below  4.0  (ca.  3.19). 

In  addition  to  chemical  modifications  of  glass,  changes  in 
the  micro  structure  of  this  material  can  be  brought  about  by 
mechanical  loading/deformation  (typically  requiring  several 
GPa  pressure  levels).  Specifically,  high  pressure  may  result  in  a 
reorganization  of  the  atomic  network  (phase  change)  in  the 
form  of  changes  to  the  coordination  of  the  network  forming 
cations.  These  phase  changes  can  be  of  first  order,  which  are 
characterized  by  the  formation  a  distinct  high-pressure  phase  at 
a  nominally  constant  pressure,  or  they  may  be  of  second  order 
which  are  phase  changes  which  involve  a  continuous  morphing 
of  the  original  phase  into  the  final  high-pressure  phase  over  a 
range  of  pressures.  These  phase  transformations  may  be 
associated  with  significant  volume  changes  and,  since  phase- 
transformation  induced  energy  absorption  is  a  well-documented 
phenomenon  responsible  for  high  toughness  levels  in  TRIP 
steels  and  partially  stabilized  crystalline  ceramics,  it  is  of 
interest  to  the  present  study. 

The  phase  transformations  of  glass  which  are  investigated  in 
the  present  work  occur  in  the  range  of  3  to  5  GPa  and  are 
thought  to  be  of  first  order.  This  pressure  range  was  chosen  for 
the  current  investigation  as  it  is  consistent  with  the  levels  seen 
in  glass  structures  subject  to  ballistic/blast  loading  and  is 
associated  with  relatively  modest  (3-7%)  volume  changes. 
These  phase  transformations  should  not  be  confused  with  those 
seen  to  take  place  at  substantially  higher  pressures  (ca. 
>20  GPa),  which  are  considered  to  be  of  second  order  and 
associated  with  substantially  larger  volume  reductions  and  with 
the  formation  of  stishovite,  an  octahedrally  coordinated  glass 
phase. 


3.  Molecular-Level  Analysis  of  Soda-Lime  Glass 

As  mentioned  earlier,  molecular-level  computational  meth¬ 
ods  have  been  employed  in  the  present  work  to  investigate 
shock  formation  and  propagation  in  soda-lime  glass.  Within 
these  methods,  all  atoms,  ions  and  bonds  are  explicitly 
accounted  for  and  molecular  mechanics,  dynamics  or  Monte 
Carlo  algorithms  are  used  to  quantify  the  behavior  of  the 
material  under  investigation. 

While  ab  initio  quantum  mechanics  methods  have  the 
advantage  over  the  molecular-level  methods  since  they  do  not 
require  parameterization,  they  have  a  serious  short-coming. 
Namely,  due  to  prohibitively  high  computational  cost,  they  can 
be  currently  employed  only  for  systems  containing  no  more 
than  a  few  hundred  atoms/particles.  As  will  be  shown  below, 
while  ab  initio  quantum-mechanics  calculations  are  not  directly 
used  in  the  present  work,  some  of  the  computational  ab  initio 
quantum  mechanics  results  are  used  in  the  parameterization  of 
the  material  model  at  the  molecular  length/time  scale.  Utility  of 
the  molecular-level  computational  results  is  greatly  dependent 
on  accuracy  and  fidelity  of  the  employed  force  field  (the 
mathematical  expression  which  describe  various  bonding  and 
non-bonding  interaction  forces  between  the  constituents  of  the 
molecular-level  model).  In  the  present  work,  the  so-called 
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“Condensed-phase  Optimized  Molecular  Potentials  for  Atom¬ 
istic  Simulation  Studies”  (COMPASS)  force  field  is  used  (Ref 
22,  23).  This  highly  accurate  force  field  is  of  an  ab  initio  type 
since  most  of  its  parameters  were  determined  by  matching  the 
predictions  made  by  the  ab  initio  quantum  mechanics  calcu¬ 
lations  to  the  condensed-matter  experimental  data.  Hence,  it 
should  be  recognized  that  the  COMPASS  force  field  is  a  prime 
example  of  how  the  highly  accurate  results  obtained  on  one 
length/time  scale  (quantum  mechanic/electronic,  in  the  present 
case)  and  the  experimental  data  can  be  combined  to  parame¬ 
terize  material  models  used  at  coarser  length/time  scale  (the 
molecular  length/time  scale,  in  the  present  case). 

Formulation  of  a  molecular-level  simulation  problem 
requires,  at  a  minimum,  specification  of  the  following  three 
aspects:  (a)  a  molecular-level  computational  model  consisting 
of  atoms,  ions,  functional  groups,  and/or  molecules;  (b)  a  set  of 
force  field  functions;  and  (c)  a  computational  method(s)  to  be 
used  in  the  simulation.  More  details  of  these  three  aspects  of  the 
molecular-level  modeling  and  simulation  of  soda-lime  glass  are 
provided  below. 


3.1  Computational  Model 

At  the  molecular  level,  soda-lime  glass  is  modeled  as  a 
discrete  material  consisting  of:  (a)  silicon  (Si)  and  oxygen  (O) 
atoms  mutually  bonded  via  a  single  covalent  bond  and  forming 
a  connected,  non-structured/amorphous  network  of  silica 
(Si044-)  tetrahedra;  (b)  oxygen  anions  (O2-)  attached  as 
terminal  functional  groups  to  the  fragmented  silica  tetrahedra 
network;  and  (c)  sodium  cations  (Na+)  dispersed  between 
fragmented  silica  tetrahedra  networks  and  ionically  bonded  to 
the  oxygen  anions. 

While  glass  is  an  amorphous  material  and  does  not  possess 
any  long-range  regularity  in  its  atomic/molecular  structure, 
modeling  of  bulk  behavior  of  glass  is  typically  done  at  the 
molecular  level  by  assuming  the  existence  of  a  larger  (amor¬ 
phous)  unit  cell.  Repetition  of  this  cell  in  the  three  orthogonal 
directions  (the  process  also  known  as  application  of  the 
“periodic  boundary  conditions”)  results  in  the  formation  of  an 
infinitely  large  bulk-type  material.  This  procedure  was  adopted 
in  the  present  work. 

The  parallelopiped-shaped  computational  cell  used  in  the 
present  work  contained  2304  particles  (Si-512,  Na+-512,  O- 
256  and  0-1024).  The  computational  cell  edge-lengths  were 
initially  set  to  a  =  7.471  nm  and  b  =  c  =  1.868  nm,  yielding  a 
soda-lime  glass  initial  nominal  density  of  2.613  g/cm3.  The 
three  edges  (a,  b ,  and  c )  of  the  cell  were  aligned,  respectively, 
with  the  three  coordinate  axes  ( x ,  y,  and  z). 

To  create  the  initial  particle  configuration  in  the  unit  cell,  the 
Visualizer  (Ref  24)  program  from  Accelrys  was  first  used  to 
construct  a  short  silica-chain  fragment.  The  fragment  was  then 
“grown”  by  a  duplicate-and-attach  process  using  the  same 
program.  The  resulting  silica  network  (along  with  additional 
sodium  cations  and  oxygen  anions)  was  next  used  within  the 
Amorphous  Cell  program  (Ref  25)  from  Accelrys  to  randomly 
populate  the  computational  cell  while  ensuring  that  the  target 
material  density  of  2.613  g/cm3  was  attained.  An  example  of  a 
typical  molecular-level  topology  within  a  single  unit  cell  is 
displayed  in  Fig.  1. 

3.2  Force  Fields 


Fig.  1  (a)  The  computational  unit  cell  for  soda-lime  glass  molecu¬ 

lar-level  simulations  used  in  the  present  work;  and  (b)  an  example  of 
the  local  atomic  structure 


which  describe,  in  an  approximate  manner,  the  potential  energy 
hyper- surface  on  which  the  atomic  nuclei  move.  In  other  words, 
the  knowledge  of  force  fields  enables  determination  of  the 
potential  energy  of  a  system  in  a  given  configuration.  In 
general,  the  potential  energy  of  a  system  of  interacting  particles 
can  be  expressed  as  a  sum  of  the  valence  (or  bond),  Ew alence, 
cross-term,  7/cr0ss-term?  and  non-bond,  7/n0n-bond?  interaction 
energies  as: 

//total  —  -^valence  T  ^cross-term  T  //non-bond  (Eq  1 ) 


The  valence  energy  generally  includes  a  bond  stretching 
term,  E^ond?  a  two-bond  angle  term,  Wangle?  a  dihedral  bond- 
torsion  term,  7/t0rsion?  an  inversion  (or  an  out-of-plane  interac¬ 
tion)  term,  Toop,  and  a  Urey-Bradlay  term  (involves  interactions 
between  two  atoms  bonded  to  a  common  atom),  EjjB,  as 


(Eq  2) 


The  cross-term  interacting  energy,  7/cross-term?  accounts  for 
the  effects  such  as  bond  length  and  angle  changes  caused  by  the 
surrounding  atoms  and  generally  includes:  stretch-stretch 
interactions  between  two  adjacent  bonds,  //bond-bond?  stretch- 
bend  interactions  between  a  two-bond  angle  and  one  of  its 
bonds,  ^bond-angle?  bend-bend  interactions  between  two  valence 
angles  associated  with  a  common  vertex  atom,  //angie-angie? 
stretch-torsion  interactions  between  a  dihedral  angle  and  one  of 
its  end  bonds,  Ee nd  bond-torsion?  stretch-torsion  interactions 
between  a  dihedral  angle  and  its  middle  bond,  EmiddiQ  bond- 
torsion?  bend-torsion  interactions  between  a  dihedral  angle  and 
one  of  its  valence  angles,  //angie-torsion?  and  bend-bend-torsion 
interactions  between  a  dihedral  angle  and  its  two  valence 
angles,  //angie-angie-torsion?  terms  as. 


Across— term  ~//bond— bond  T~  Wangle— angle  T  //bond— angle 
T  //end_bond— torsion  T  -^middle-bond— torsion 
T  -Wangle— torsion  T  //angle— angle— torsion 


(Eq  3) 


As  stated  above,  the  behavior  of  a  material  system  at  the 
molecular-level  is  governed  by  the  appropriate  force  fields 


The  non-bond  interaction  term,  7/n0n-bond?  accounts  for  the 
interactions  between  non-bonded  atoms  and  includes  the  van 
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der  Waals  energy,  Ev dw,  and  the  Coulomb  electrostatic  energy, 

^Coulomb? 

^non-bond  —  -^vdW  T~  ^Coulomb  (Eq  4) 

As  mentioned  earlier,  the  present  molecular-level  analysis  of 
soda-lime  glass  employs  the  COMPASS  (Ref  22,  23)  force  field 
for  various  bond  and  non-bond  interaction  energies  appearing 
in  Eqs  1-4.  A  summary  of  the  COMPASS  force-field  functions 
can  be  found  in  our  previous  work  (Ref  26). 

3.3  Computational  Method 

Both  molecular  statics  and  molecular  dynamics  simulations 
were  employed  in  the  present  work.  Within  the  molecular 
statics  approach,  the  unit-cell  potential  energy  (as  defined  by 
Eqs  1-4)  is  minimized  with  respect  to  the  position  of  the 
constituent  particles/atoms.  The  potential  energy  minimization 
within  Discover  (Ref  27)  (the  atomic  simulation  program  from 
Accelrys  used  in  the  present  work)  is  carried  out  by  combining 
the  Steepest  Descent,  Conjugate  Gradient,  and  the  Newton’s 
minimization  algorithms.  These  algorithms  are  automatically 
inactivated/activated  as  the  atomic  configuration  is  approaching 
its  energy  minimum  (i.e.,  the  Steepest  Descent  method  is 
activated  at  the  beginning  of  the  energy  minimization  procedure 
while  the  Newton’s  method  is  utilized  in  the  final  stages  of  this 
procedure). 

As  will  be  discussed  in  greater  detail  in  Sect.  4,  molecular 
statics  is  employed  to  determine  the  state  of  the  material  swept 
by  a  shock.  As  will  be  shown,  this  procedure  is  based  on  the 
use  of  (bonding  and  non-bonding)  potential  energy  components 
and  neglects  shock-induced  changes  in  the  (configurational) 
entropy  of  the  system.  To  assess  the  consequence  of  this 
simplification,  the  approach  described  in  Ref  28  was  consid¬ 
ered.  This  approach  defines  a  dimensionless  parameter  and 
states  that  when  this  parameter  is  significantly  smaller  than 
unity,  entropy  effects  can  be  neglected.  Unfortunately,  detailed 
temperature  and  pressure  dependencies  of  the  material  mechan¬ 
ical  response  of  soda-lime  glass  needed  to  evaluate  this 
parameter  were  not  available  with  sufficient  fidelity.  Hence, 
the  results  obtained  by  the  application  of  this  procedure,  which 
suggest  that  the  entropy  effects  may  not  be  highly  critical, 
cannot  be  accepted  with  a  high  level  of  confidence. 

Within  the  molecular  dynamics  approach,  gradient  of  the 
potential  energy  with  respect  to  the  particle  positions  is  first 
used  to  generate  forces  acting  on  the  particles  and,  then,  the 
associated  Newton’s  equations  of  motion  (for  all  particles)  are 
integrated  numerically  to  track  the  temporal  evolution  of  the 
particle  positions.  Both  the  equilibrium  and  non-equilibrium 
molecular  dynamics  methods  are  employed  in  the  present  work. 
Within  the  equilibrium  molecular-dynamics  methods,  the 
system  under  consideration  is  coupled  to  an  (external)  envi¬ 
ronment  (e.g.,  constant  pressure  piston,  constant  temperature 
reservoir,  etc.)  which  ensures  that  the  system  remains  in 
equilibrium  (i.e.,  the  system  is  not  subjected  to  any  thermody¬ 
namic  fluxes).  As  will  be  discussed  in  next  section,  NVT  (where 
N  is  the  (fixed)  number  of  atoms,  V,  the  computational  cell 
volume  (also  fixed),  and  T  (298  K)  is  the  temperature) 
equilibrium  molecular  dynamics  is  employed  in  the  first  stage 
of  the  shock  generation  procedure.  In  addition,  as  will  be 
discussed  in  Sect.  4,  NVE  ( E  is  the  total  energy)  equilibrium 
molecular  dynamics  is  also  employed  during  determination  of 
the  shock  Hugoniot.  Within  non-equilibrium  molecular  dynam¬ 
ics,  the  system  is  subjected  to  large  perturbations  (finite 


changes  in  the  axial  parameter  of  the  computational  cell,  in  the 
present  case)  which  create  a  thermodynamic  flux  (i.e.,  the  flux 
of  energy  and  momentum,  in  the  present  case).  More  details 
regarding  the  use  of  Discover  to  carry  out  molecular  statics  and 
molecular  dynamics  analyses  can  be  found  in  our  prior  work 
(Ref  6). 

3.4  Shock- Wave  Generation 

To  generate  a  planar  shock  (or  more  precisely  a  pair  of 
planar  shocks)  within  the  computational  cell,  the  following 
procedure  is  employed: 

(a)  At  the  beginning  of  the  analysis,  a  “sufficiently  long” 
NVT  molecular  dynamics  simulation  is  carried  out  to 
equilibrate  the  system/material. 

(b)  The  shock  is  then  initiated  (and  driven)  by  continuously 
contracting  the  computational  cell  x-direction  lattice 
parameter  a  as: 

a(t)  —  a(t  =  0)  —  2upt  (Eq  5) 

where  t  denotes  time,  uv  is  the  so-called  “piston”  velocity  (or 
equivalently  the  particles  upstream  velocity)  in  the  x-direc- 
tion.  uv  is  varied  over  a  range  between  187.5  and  1500  m/s 
to  simulate  the  generation  and  propagation  of  shock  of  vari¬ 
ous  strengths.  Meanwhile,  computational-cell  transverse  lat¬ 
tice  parameters  b  and  c  are  kept  constant  to  obtain  planar 
(uniaxial-strain)  shock  conditions.  In  this  process,  the  compu¬ 
tational  cell  faces  normal  to  the  shock  propagation-direction 
behave  very  similarly  to  the  impact  surface  of  a  plate-like  tar¬ 
get  subjected  to  a  so-called  symmetric  “flyer-plate”  impact 
test  (Ref  29).  The  procedure  employed  here  generates  a  pair 
of  shock  waves  which  propagate,  at  a  shock  speed  US9  from 
the  cell  boundaries  toward  its  center.  As  schematically  shown 
in  Fig.  2,  these  shock  waves  leave  behind  a  “shocked”  mate¬ 
rial  state  characterized  by  a  higher  material  density  (as  well 
as  internal  energy,  temperature,  stress,  particle  velocity,  and 
entropy). 


Fig.  2  A  schematic  of  the  generation  of  a  pair  of  shocks  in  a 
molecular-level  system  via  the  process  of  computational-cell  parame¬ 
ter  contraction 
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The  aforementioned  procedure  for  shock-wave  generation 
and  the  subsequent  molecular  statics/dynamics  analyses  are 
carried  out  through  the  use  of  a  Discover  input  file  (Ref  27) 
which  is  written  in  a  Basic  Tool  Command  Language  (BTCL). 
This  enabled  the  use  of  a  scripting  engine  that  provides  very 
precise  control  of  simulation  jobs,  e.g.,  a  cell  deformation  to  be 
carried  out  in  small  steps  each  followed  by  a  combined  energy- 
minimization/molecular-dynamics  simulation  run. 

3.5  Problem  Formulation 

The  problem  addressed  in  the  present  work  involved 
generation  of  shock  waves  of  different  strengths  (using  the 
aforementioned  computational  cell  parameter  contraction  meth¬ 
od),  determination  of  the  associated  shock  Hugoniot  relations 
and  identification  and  elucidation  of  the  main  molecular-level 
inelastic-deformation/energy  dissipation  processes  taking  place 
at  or  in  the  vicinity  of  the  shock  front.  The  procedure  for  shock- 
wave  generation  was  presented  in  the  previous  section. 

As  far  as  the  shock  Hugoniot  determination  is  concerned,  it 
entailed  the  knowledge  of  the  shock-wave  profiles  (and  their 
temporal  evolution)  for  the  axial  stress,  material  density, 
particle  velocity,  internal  energy,  and  temperature.  The  latter 
are  obtained  by  lumping  particles/atoms  and  their  (bond  and 
non-bond)  potential  and  kinetic  energy  contribution,  into  fixed- 
width  bins,  in  the  order  of  their  axial  coordinates.  As  will  be 
shown  in  the  next  section,  two  types  of  bins  are  used:  (a)  a 
Lagrangian-type  which  is  fixed  to  the  initial/reference  state  of 
the  computational  cell  and  (b)  a  moving-type  which  is  attached 
to  the  advancing  shock  front. 

Identification  of  the  molecular-level  inelastic-deformation/ 
energy  dissipation  processes  entailed  a  close  examination  of  the 
changes  in  a  material  bond  structure  and  topology  caused  by 
the  passage  of  the  shock. 


4.  Results  and  Discussion 

4.1  Shock- Wave  Observation  and  Structure 
Characterization 

An  example  of  the  typical  results,  obtained  in  the  present 
work,  pertaining  to  material  molecular-level  microstructure/ 
topology  evolution  caused  by  a  continuous  axial  contraction  of 
the  computational  cell  is  displayed  in  Fig.  3(a)-(d).  The  results 
displayed  in  this  figure  clearly  reveal  the  generation  of  a  pair  of 
planar  shock  waves  at  the  two  opposing  y-z  faces  of  the 
computational  cell,  Fig.  3(a),  and  their  subsequent  propagation 
toward  the  center  of  the  cell,  Fig.  3(b)-(d).  An  approximate 
location  of  the  center  point  of  the  two  shocks  is  indicated  using 
arrowheads  in  Fig.  3(a)-(d).  The  results  displayed  in  Fig.  3(a)- 
(d)  show  that  the  shock  waves  remain  fairly  planar  during  their 
motion.  While  the  two  shocks  collide  at  a  later  simulation  time, 
a  shock-wave  interaction  investigation  is  beyond  the  scope  of 
the  present  work. 

While  the  results  displayed  in  Fig.  3(a)-(d)  provide  clear 
evidence  for  the  formation  and  propagation  of  a  pair  of 
opposing  planar  shock  waves,  they  do  not  offer  any  information 
about  the  structure/shape  of  the  shock-wave  front  or  about  the 
state  of  the  (upstream)  material  swept  by  the  shock.  The  latter 
aspects  of  shock-wave  generation  and  propagation  within  soda- 
lime  glass  are  addressed  in  the  remainder  of  this  section  and  in 
the  subsequent  sections. 


Fig.  3  Temporal  evolution  of  the  molecular  level  material  micro¬ 
structure  accompanying  generation  and  propagation  of  a  pair  of  pla¬ 
nar  shocks  in  soda-lime  glass 


To  reveal  the  structure/shape  of  the  shock  wave,  the  method 
of  (Lagrangian)  bins  described  in  Sect.  3.4  is  employed.  In  this 
case,  the  bins  are  fixed  in  the  (initial)  reference  configuration  of 
the  computational  cell.  In  other  words,  the  same  atoms  are 
associated  with  a  given  bin  throughout  the  entire  molecular 
dynamics  simulation.  Examples  of  the  typical  results  obtained 
through  the  use  of  this  method  are  displayed  in  Fig.  4(a)-(b). 
The  results  displayed  in  this  figure  are  obtained  under  identical 
conditions  except  for  the  rate  of  axial  contraction  of  the 
computational  cell  (100%  higher  in  the  case  of  Fig.  4(b)).  In 
these  figures,  particle  velocities  at  different  simulation  (i.e., 
post-shock  wave  generation)  times  are  plotted  against  the 
Lagrangian  bin  center  x-location  (for  the  shock  propagating  to 
the  right,  only).  Brief  examination  of  the  results  displayed  in 
Fig.  4(a)-(b)  reveals  that: 

(a)  two  shock  waves  are  generated  (only  the  right-propagat¬ 
ing  shock  is  shown  though)  at  the  computational  cell 
faces  normal  to  the  x-direction.  These  shocks  subse¬ 
quently  propagate  towards  the  computational-cell  center; 

(b)  after  a  brief  transient  period,  the  shocks  appear  to  attain 
steady  wave  profiles  (i.e.,  a  time-invariant  profile  within 
a  reference  frame  which  is  attached  to,  and  moves  with, 
the  shock  front);  and 

(c)  both  the  particle  velocity  and  the  shock  speed  increase 
with  the  computational-cell  contraction  rate.  It  should  be 
noted  that  the  curves  bearing  the  same  numerical  label 
in  Fig.  4(a)-(b)  correspond  to  the  same  simulation  time. 

It  should  also  be  noted  that  no  thermostat  was  used  in  the 
present  non-equilibrium  molecular  dynamics  simulations,  so 
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(b)  Distance  in  x-direction,  nm 


Fig.  4  Temporal  evolution  of  the  particle  velocity  associated  with 
the  propagation  of  two  approaching  shock  waves.  The  results  in  (b) 
are  associated  with  a  50%  higher  computational-cell  axial  contrac¬ 
tion  rate 


that  the  steady-wave  shock  profile  is  a  natural  consequence  of  a 
balance  between  the  continuous  supply  of  momentum  to  the 
system  (through  the  continuous  computational  cell  axial 
contraction)  and  the  observed  lateral  motion  of  the  atoms  in 
the  continuously  enlarged  upstream  material  domain  swept  by 
the  shock.  In  general,  the  use  of  a  thermostat  modifies  the 
(F  =  ma  Newtonian- type)  equations  of  motion  solved  within 
the  molecular  dynamics  simulations  by  the  introduction  of  a 
(velocity-proportional)  viscous-dissipation  term.  It  is  well 
established  that,  when  shock  formation  and  propagation  is 
analyzed  within  a  continuum  framework,  the  use  of  a  viscous- 
dissipation  term  is  mandatory  for  the  attainment  of  a  steady- 
wave  shock  profile.  This  fact  has  often  been  used  as  a 
justification  for  the  use  of  a  (local  or  global)  thermostat  within 
molecular-level  simulations  of  shock-wave  formation/propaga- 


Fig.  5  The  Hugoniot  relation  pertaining  to  the  particle  velocity 
dependence  of  the  shock  speed 

tion.  While  such  practices  greatly  facilitate  the  attainment  of  a 
steady  shock,  they  cannot  be  readily  defended  since  shock 
formation  and  propagation  is  generally  considered  to  be  an 
adiabatic  (no  system/surrounding  energy  exchange)  process  due 
to  the  near-instantaneous  material  transition  to  the  shocked 
state.  In  addition,  shock  generation/propagation  is  a  non- 
isentropic  process  due  to  the  attendance  of  various  energy 
dissipation  mechanisms.  These  were  the  determining  factors  in 
the  decision  to  not  use  a  thermostat  in  the  present  work.  Due  to 
the  lack  of  a  thermostat,  the  molecular  dynamics  simulations 
employed  in  the  present  work  can  be  characterized  as  being  of  a 
non-equilibrium  type.  It  is  interesting  to  point  out  that  despite 
the  fact  that  no  viscous-dissipation  term  was  added  to  the 
Newton’s  equations  of  motion,  the  results  displayed  in 
Fig.  4(a)-(b)  show  some  of  the  defining  features  of  shock 
waves  when  they  are  analyzed  in  the  continuum-level  simula¬ 
tions  (in  the  presence  of  viscous  dissipation).  Specifically,  a 
steady  shock  is  generated  and  the  shock  width  decreases  with 
an  increase  in  the  shock  strength. 

4.2  Determination  of  Shock  Hugoniot  Relations 

The  results  presented  in  Fig.  4(a)-(b)  reveal  the  steady- 
shock  profile  and  can  be  used  to  find  a  functional  relation 
between  the  shock  speed,  Us,  and  the  particle  velocity,  up.  This 
was  done  in  the  present  work  and  the  results  obtained  are 
displayed  in  Fig.  5.  The  Us  vs.  up  relation  is  one  of  the  so- 
called  shock  Hugoniot  relations.  In  general,  the  Hugoniot  can 
be  defined  as  a  locus  of  shocked-material  states  in  a  stress/ 
pressure,  energy  density,  mass  density  (or  specific  volume), 
temperature,  particle  velocity,  and  shock  speed  space  which  are 
associated  with  (or  “centered  at”)  a  given  initial/reference 
material  state.  The  Us  vs.  up  relation  mentioned  above  is  a 
simple  projection  of  the  Hugoniot  to  the  Us  —  uv  plane.  In  the 
case  of  planar  shocks,  of  interest  in  the  present  work,  the  other 
commonly  used  Hugoniot  relations  include:  axial  stress,  t\ i;  vs. 
density,  p  (or  specific  volume,  v  =  1/p);  (mass-based)  internal 
energy  density,  e,  vs.  p  (or  v);  tn  vs.  up  and  temperature  T  vs.  p 
(or  v).  These  relations  were  determined  in  the  present  work 
using  two  distinct  methods: 
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(a)  The  first  method  is  based  on  the  three  so-called  “jump 
equations”  which  are  defined  as: 

p-Us=  p+(Us-up)  (Eq  6) 

tu  +  p- U2S  =  t+  +  p+  (Us  -  «p)2  (Eq  7) 

+  p  +  0.5 U2S  =  e+  +  +  0.5  (Us  -  up)Z  (Eq  8) 

These  equations  relate  the  known  downstream  material  states 
(denoted  by  a  superscript  “  —  ”)  and  the  unknown  upstream 
material  states  (denoted  by  a  superscript  “+”)  associated  with 
the  shock  of  a  given  strength  (as  quantified  by  the  shock 
speed  or  the  downstream-to-upstream  particle  velocity  jump). 
These  equations  are  next  combined  with  the  previously  deter¬ 
mined  Us  vs.  up  relation  and  the  prescribed  (shock- strength 
defining  quantity)  Us  or  up  to  solve  for  the  unknown  up¬ 
stream  material  states.  It  should  be  noted  that  this  method  en¬ 
ables  determination  of  only  material  mechanical  state 
variables  (tn,  e ,  v  ( =l/p ),  US9  and  up).  To  obtain  temperature, 
a  separate  set  of  equilibrium  NVE  (E — total  energy  of  the 
system)  molecular  dynamics  simulations  is  carried  out.  In 
each  case,  a  local  computational  sub-cell  is  defined  containing 
only  the  upstream  (shocked)  material.  The  number  of  parti¬ 
cles,  the  volume  of  the  sub-cell  and  its  total  internal  energy 
are  all  maintained  constant.  The  associated  “equilibrium” 
temperature  is  then  calculated  using  the  time-averages  of  the 
atomic  velocities  (see  Eq  10);  and 


(b)  Time  averages  of  the  atomic  positions,  rh  velocities,  vh 
and  interaction  forces,  f  (i  is  the  atomic  label)  are  used 
to  compute  the  unknown,  local  (bin-based)  thermo¬ 
mechanical  quantities  using  the  following  standard  sta¬ 
tistical  thermodynamic  relations: 


1 

Vbm 


Nhin 

Em‘ 

i= 1 


avg 


T  = 
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3A^bin&b 
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«bi„ 
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tu  =  -pr—  (  Vbin kbT  +  rt  ®fi 


(Eq  9) 


(Eq  10) 


(Eq  11) 


E  =  (  E’Total  +  ^  \  EmiVi  ’  Vi 


(Eq  12) 


where  subscript  “avg”  denotes  time  averaging,  Nhin  and  Vbin 
the  number  of  atoms  within  and  the  volume  of  the  bin, 
respectively,  kh  is  the  Boltzmann’s  constant,  ETota\  is  given  by 
Eq  1,  while  “•”  and  “  ®”  indicate  dot  product  and  tensorial 
product  operators,  respectively. 

It  should  be  noted  that,  in  this  case,  the  bins  were  defined 
within  a  reference  frame  which  is  attached  to  (and  moves  with) 
the  steady-shock  front.  Clearly,  in  this  case  different  atoms  may 
reside  within  a  given  bin  at  different  simulation  times.  On  the 
other  hand,  the  bins  correctly  collect  the  information  about  the 
atoms  (temporarily)  residing  in  a  given  segment  of  the  steady- 
shock  profile.  In  other  words,  time  averages  are  calculated  not 


for  a  fixed  assembly  of  atoms,  but  rather  for  a  transient  set  of 
atoms  associated  with  a  given  moving  bin.  It  should  be  also 
noted  that  since  one  of  the  main  objectives  of  the  present  work 
was  determination  of  the  Hugoniot  relations,  only  the  data 
pertaining  to  the  bins  located  in  the  fully  shocked  upstream 
region  are  collected  and  analyzed  (for  different  shock-strength 
conditions).  In  other  words,  the  data  collected  by  the  bins 
located  within  the  shock  profile  and  downstream  of  the  shock 
are  ignored. 

The  results  of  these  two  procedures  are  displayed  in 
Fig.  6(a)-(d).  In  this  figure,  two  cases  are  shown  and  labeled 
as  “Method  (a)”  or  “Method  (b)”  in  order  to  indicate  the 
method  used  for  generation  of  the  corresponding  results.  It  is 
apparent  that  the  two  methods  yield  consistent  results. 

The  Hugoniot  relations  displayed  in  Fig.  5  and  6(a)-(d)  are 
typically  used  within  a  continuum-level  computational  analysis 
of  shock-wave  generation/  propagation  in  two  ways: 

(a)  They  are  directly  used  in  the  analysis  of  shock-wave 
propagation  under  uniaxial  strain  conditions  (Ref  30, 
31);  and 

(b)  Alternatively,  they  can  be  used  to  derive  a  continuum- 
level  material  model  which  is  consistent  with  the  mate¬ 
rial  mechanical  response  under  high-rate,  large-strain, 
high-pressure  conditions.  Such  a  model  is  subsequently 
used  in  general  three-dimensional,  non-linear  dynamics 
computational  analyses  (Ref  32). 

In  our  recent  work  (Ref  6),  it  was  shown  that  the  Hugoniot 
relations  can  be  generated  by  converting  the  corresponding 
isotherms  (obtained  via  quasi-static,  molecular-level  mechan¬ 
ical  tests).  This  procedure  was  found  to  be  associated  with  a 
number  of  challenges  (e.g.,  a  particular  form  of  the  equation  of 
state  had  to  be  assumed,  several  material  properties/relations 
had  to  be  assessed  independently,  etc.).  Most  of  these 
challenges  were  not  encountered  in  the  present  work  since 
the  Hugoniot  relations  are  derived  more  directly  from  the 
molecular-level  computational  results. 

4.3  Shock-Induced  Material-State  Changes 

The  results  presented  and  discussed  in  the  previous  sections 
clearly  revealed  the  formation  and  propagation  of  planar  shocks 
in  soda-lime  glass  and  enabled  formulation  of  the  appropriate 
shock  Hugoniot  relations.  In  the  present  section,  a  more 
detailed  investigation  is  carried  out  of  the  molecular-level 
material  microstructure  in  the  wake  of  a  propagating  planar 
shock. 

It  should  be  first  realized  that  the  analysis  of  material 
micro  structure  and  its  evolution  due  to  shock  loading  in  soda- 
lime  glass  is  quite  challenging  due  to  the  absence  of  a  crystal 
structure  in  this  material.  Namely,  when  molecular-level 
simulations  of  shock  generation/propagation  are  carried  out  in 
(nearly  perfect)  single  crystal  solids  (Ref  1 7,  20),  the  presence 
of  a  crystal  lattice  greatly  facilitates  the  investigation  of 
deviations  from  the  long  range  order  (i.e.,  formation  of  various 
point,  line  and  planar  defects)  and  the  nature  of  the  associated 
inelastic  deformation  processes.  Soda- lime  glass  is,  on  the  other 
hand,  an  amorphous  material  in  its  initial  condition  and  remains 
so  after  being  subjected  to  shock  loading.  To  address  the 
challenge  of  material  micro  structure  characterization  and  its 
changes  resulting  from  shock  loading,  two  micro  structural 
parameters  are  monitored  in  the  present  work:  (a)  the  extent  of 
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Fig.  6  (a)  Axial  stress  vs.  specific  volume;  (b)  energy  vs.  specific  volume;  (c)  axial  stress  vs.  particle  velocity;  and  (d)  temperature  vs.  specific 

volume  Hugoniot  relations.  Please  see  text  for  explanation  of  the  “Method  (a)”  and  “Method  (b)” 


lateral  motion  of  the  atoms  at  the  shock  front;  and  (b)  the  extent 
of  changes  to  the  Si  and  O  atomic  coordination  (i.e.,  bonding 
structure). 

In  the  case  of  single-crystalline  solids,  previous  molecular- 
level  shock-wave  formation/calculation  work  (Ref  17,  20) 
established  that  the  steady-wave  condition  is  attained  not  as  a 
result  of  (velocity-dependent)  viscous  dissipation  (as  is  the  case 
for  shocks  in  fluids),  but  rather  as  a  result  of  lateral  atomic 
motions  which  result  in  inelastic  deformation  (permanent 
slippage  of  crystal  planes  and  the  formation  of  crystal  defects). 
Similar  lateral  motions  of  the  atoms  at,  or  in  the  vicinity  of,  the 
shock  front  are  also  observed  in  the  present  work.  Since  these 
motions  result  in  a  permanent  change  in  atomic  nearest 
neighbor  coordination  (as  well  as  in  a  permanent  volume 
change  of  the  shocked  material),  they  can  be  considered  as 
inelastic  strain  producing.  However,  as  stated  earlier,  the 
amorphous  nature  of  soda-lime  glass  precludes  a  more  detailed/ 
quantitative  description  of  the  nature  of  the  atomic  lateral 


motion-induced  inelastic-deformation  process(s)  or  the  associ¬ 
ated  micro  structural  defect  formation.  On  the  other  hand,  the 
accompanying  changes  in  the  atomic  coordination  and  the  bond 
structure  can  be  readily  identified  (as  will  be  demonstrated 
below). 

Examples  of  the  shock-loading  induced  changes  in  the 
atomic  coordination  and  bond  structure  of  soda-lime  glass  are 
given  in  Fig.  7(a)-(d)  and  8(a)-(b).  To  aid  in  visualization/ 
interpretation  of  the  microstructural/topological  changes  expe¬ 
rienced  by  soda- lime  glass  during  shock  loading,  only  a  20-30 
atom  exemplary  region  of  the  upstream  material  is  analyzed  in 
these  figures.  In  addition,  sodium  cations  are  not  displayed. 

In  Fig.  7(a)  and  (c),  x-y  and  y-z  projections  are  given  of  the 
material  region  under  consideration  in  its  initial  state.  The 
corresponding  projections  for  the  same  material  region  after 
the  passing  of  a  shock  are  displayed  in  Fig.  7(b)  and  (d).  A 
closer  examination  of  the  microstructure/topology  results  dis¬ 
played  in  Fig.  7(a)-(d)  provides  clear  evidence  for  the  formation 
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Fig.  7  Development  of  five-fold  coordinated  silicon  atoms  in  soda- 
lime  glass  under  shock  loading:  (a)  material  initial  state,  x-y  projec¬ 
tion;  (b)  post-shock  state  of  the  same  material  region  as  in  (a);  (c) 
and  (d)  correspond  respectively  to  the  y-z  projection  of  (a)  and  (b) 

of  five-fold  coordinated  silicon  atoms  (circled  for  improved 
clarity).  In  addition,  a  comparison  of  the  results  displayed  in 
Fig.  7(a)-(b)  reveals  the  extent  (partly  irreversible)  compaction 
caused  by  shock  loading  of  soda-lime  glass. 

Additional  typical  changes  in  the  soda-lime  glass  molecular- 
level  microstructure  induced  by  shock  loading  are  displayed  in 
Fig.  8(a)-(b).  The  results  displayed  in  Fig.  8(a)  pertain  to  the 
material  initial  state  while  those  shown  in  Fig.  8(b)  correspond 
to  the  shocked  state  of  the  same  material  region.  A  close 
examination  of  the  results  displayed  in  this  figure  reveals  the 
formation  of  smaller  Si-0  rings.  For  improved  clarity,  these 
rings  are  highlighted  in  purple  and  yellow  in  Fig.  8(b). 

It  should  be  noted  that  the  molecular-level  microstructural 
changes  described  above  are  a  manifestation  of  the  transition  of 
soda-lime  glass  to  a  state  that  is  energetically  favored  at  high 
shock- induced  stresses  (large  densities).  This  finding  is  con¬ 
sistent  with  the  experimental  observation  of  Alexander  et  al. 
(Ref  32),  who  reported  the  formation  of  stishovite-like  domains 
containing  six-fold  coordinated  silicon  atoms  in  soda-lime  glass 
at  ca.  30  GPa  shock-induced  axial  stress  levels. 

It  should  be  also  recognized  that  shock  loading  leads  to  a 
permanent  2-4%  (shock  strength-dependent)  increase  in  the 
material  density.  This  material  densification  process  along  with 
the  aforementioned  processes  leading  to  an  increase  in  the 
silicon-atom  coordination  and  smaller  Si-0  ring  formation  are 
all  associated  with  energy  absorption/dissipation  and,  hence, 
are  expected  to  play  an  important  role  in  the  blast/ballistic 
impact  mitigation  potential  of  soda-lime  glass. 

4.4  Final  Remarks 

Within  the  present  work,  molecular-level  computational 
methods  are  employed  to  study  various  phenomena  accompa¬ 


w 
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Fig.  8  Formation  of  smaller  Si-0  rings  in  soda-lime  glass  under 
shock  loading:  (a)  material  initial  state,  y-z  projection;  and  (b)  post¬ 
shock  state  of  the  same  material  region  as  in  (a).  Newly  formed  Si-0 
ring  structures  are  highlighted  for  clarity 

nying  shock-wave  generation  and  propagation  in  soda-lime 
glass.  It  should  be  noted  that  the  present  work  does  not  suggest 
that  molecular-level  analyses  of  shock  generation  and  propa¬ 
gation  should  replace  the  corresponding  continuum-level 
(hydrodynamic)  analyses.  The  latter  are  far  better  suited  (and 
feasible)  for  studying  the  behavior  of  real-life  engineering 
systems  (e.g.,  vehicle  windshield  subjected  to  blast  impact). 
Rather,  the  present  approach  is  highly  beneficial  relative  to  the 
identification  and  characterization  of  the  nanoscale  phenomena/ 
processes  taking  place  at  the  shock  front.  Once  these  phenom¬ 
ena/processes  are  well  understood  and  characterized  (a  formi¬ 
dable  task)  the  knowledge  gained  can  be  used  to  formulate  (and 
parameterize)  more  physically  based  material  models  suitable 
for  use  in  continuum-level  analyses. 


5.  Summary  and  Conclusions 

Based  on  the  results  obtained  in  the  present  work,  the 
following  summary  remarks  and  main  conclusions  can  be  drawn: 

1 .  Various  phenomena  accompanying  the  formation  and 
propagation  of  a  planar  shock  wave  within  soda-lime 
glass,  a  material  commonly  used  in  transparent  armor 
applications,  are  investigated  using  molecular-level  com¬ 
putational  methods. 
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2.  The  results  obtained  show  that  even  without  the  use  of  a 
viscous-dissipation-based  thermostat,  a  steady-wave  pla¬ 
nar  shock  profile  can  readily  be  established  in  this  mate¬ 
rial. 

3.  The  time-averaged  results  pertaining  to  the  atomic  posi¬ 
tions,  velocities  and  interaction  forces  are  used  to  con¬ 
struct  the  appropriate  shock  Hugoniot  relations,  the 
relations  which  define  the  locus  of  stress,  energy  density, 
mass  density,  temperature  and  particle  velocity  in  soda- 
lime  glass  swept  by  a  shock  propagating  at  a  given 
speed. 

4.  Detailed  examination  of  the  molecular-level  microstruc¬ 
ture  evolution  in  the  shock-wave  wake  is  carried  out  to 
identify  the  nature  of  energy-absorbing  and  shock  wave 
spreading  mechanisms.  The  results  revealed  that  shock 
loading  causes  extensive  changes  in  atomic  coordination 
and  the  bond  structure  as  well  as  a  2-4%  (shock  strength- 
dependent)  density  increase.  Specifically,  the  atomic  coor¬ 
dination  of  many  silicon  atoms  has  been  found  to 
increase  from  four  to  five  and  numerous  smaller  Si-0 
rings  are  observed.  These  processes  are  associated  with 
substantial  energy  absorption  and  dissipation  and  are  be¬ 
lieved  to  greatly  influence  the  blast/ballistic  impact  miti¬ 
gation  potential  of  soda-lime  glass. 
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